Imports
library(evd);
library(evdbayes);
library(coda);
library(ismev);
Lade nötiges Paket: mgcv
Lade nötiges Paket: nlme
This is mgcv 1.8-24. For overview type 'help("mgcv-package")'.
library(xts);
Lade nötiges Paket: zoo
Attache Paket: ‘zoo’
The following objects are masked from ‘package:base’:
as.Date, as.Date.numeric
library(ggplot2);
Loading the Data
load("../data/CAPE_Minder_Rychener_Malsot.RData")
load("../data/NINO34.RData")
load("../data/SRH_Minder_Rychener_Malsot.RData")
Generate PROD
prod = sqrt(cape)*srh
** Create Time Series Objects **
dates <- seq.Date(as.Date("1979-1-1"), as.Date("2015-12-31"), by=1)
feb29ix <- format(as.Date(dates), "%m-%d") == "02-29"
dates <- dates[!feb29ix]
prod_ts = xts(prod, order.by = rep(dates, each=8))
Plot of Time Series
ts_plot = autoplot(prod_ts) + xlab("Time") + ylab("PROD")
ggsave("../plots/full_time_series.pdf", plot=ts_plot, width = 6, height = 3)
ggsave("../plots/full_time_series.png", plot=ts_plot, width = 3, height = 1.5)
prod_full_monthly_max = apply.monthly(prod_ts, max)
max_ts_plot = autoplot(prod_full_monthly_max) +
xlab("Time") + ylab("PROD") +
geom_line(size=.5)
ggsave("../plots/monthly_max_series.pdf", plot=max_ts_plot, width = 6, height = 3)
ggsave("../plots/monthly_max_series.png", plot=max_ts_plot, width = 3, height = 1.5)
nino34_plot = ggplot() + geom_line(aes(index(prod_full_monthly_max), nino34)) +
xlab("Time") + ylab("NINO 3.4")
ggsave("../plots/nino34_plot.pdf", plot=nino34_plot, width = 6, height = 3)
ggsave("../plots/nino34_plot.png", plot=nino34_plot, width = 3, height = 1.5)
Group Monthly
month_names = c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
get_monthly = function(x) {
output = list()
len = nrow(x)
# Get month of element
month = time(x)
month = gsub("....-", "", month)
month = gsub("-..", "", month)
monthlist = unique(month)
for (i in 1:12) {
output[[i]] = x[month == monthlist[i],]
}
names(output) = month_names
return(output)
}
monthly_max = get_monthly(apply.monthly(prod_ts, max))
r = 2
r_monthly_max = get_monthly(apply.monthly(prod_ts, function(x) order(x, decreasing=T)[1:r]))
Plot Monthly Time Series
for (i in 1:12) {
tmp = autoplot(monthly_max[[i]]) + xlab("Time") + ylab("PROD")
ggsave(paste("../plots/monthly_max_series/", sprintf("%02d", i), "monthly_max_series.pdf", sep=""),
plot=tmp, width = 6, height = 3)
}
# get the monthly maxima of prod
m1 = as.data.frame(apply.monthly(prod_ts, max))$V1;
# produce the gumbel qq plot
gumbel_qq = function(x, title="") qqplot(qgumbel(c(1:length(x))/(length(x)+1)),
x,
main = paste("Gumbell Q-Q Plot", title),
xlab = "Theoretical Quantiles",
ylab = "Sample Quantiles") ;
gumbel_qq(m1)

#qqplot(qgumbel(c(1:length(m1))/(length(m1)+1)),m1,main = "Gumbell Q-Q Plot",xlab = "Theoretical Quantiles", ylab = "Sample Quantiles") ;
The qq plot doesn’t fit very well, especially in the lower tail. This is likely due to seasonal dependence.
Fitting GEV to the entire data
# fit gevd with MLE and produce diagnostic plots
fitmax.MLE<-fgev(m1,method="Nelder-Mead")
par(mfrow=c(2,2))
fitmax.MLE
Call: fgev(x = m1, method = "Nelder-Mead")
Deviance: 9272.599
Estimates
loc scale shape
7.519e+03 7.013e+03 4.757e-03
Standard Errors
loc scale shape
421.90201 331.43749 0.05347
Optimization Information
Convergence: successful
Function Evaluations: 171
plot(fitmax.MLE)

Poor fit, probably because the distribution isn’t stationary. This is especially visible in the Probability plot, in which the confidence band is surpassed, indicating a poor fit.
# fit gevd with Bayesian Techniques
# use the previous outputs (rounded) as initial values (use a different shape)
init<-c(1.6e4,4e3,0.1)
# arbitrary priors
mat <- diag(c(10000,10000,100))
pn <- prior.norm(mean=c(0,0,0),cov=mat)
# proposal standard deviation (found by trial and error to get 20%<acceptance rate<40%)
psd<-c(500,0.03,0.02)
# draw 3k samples from markov chain
nit = 30000
thinning = 300
post<-posterior(nit, init=init,prior=pn,lh="gev",data=m1,psd=psd)
# diagnostic plots
MCMC<-mcmc(post[1:nit %% thinning == 0, ], thin=300)
plot(MCMC)

attr(mcmc(post),"ar")
mu sigma xi total
acc.rates 0.24 0.71 0.70 0.55
ext.rates 0.00 0.00 0.01 0.00
#MCMC_stab <- mcmc(post, thin=50, start=1000)
acf(mcmc(post[1:nit %% thinning == 0, ]))

We observe that there seem to be no substantial issues from the autocorrelation plots.
apply(mcmc(post[1:nit %% thinning == 0, ]),2,mean)
mu sigma xi
236.0524535 11181.3806761 -0.1897774
apply(mcmc(post[1:nit %% thinning == 0, ]),2,sd)
mu sigma xi
98.60894002 624.46021143 0.03409059
** Fit with MLE for months separately**
#monthly_fits = lapply(monthly_max,
# function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
monthly_fits = list()
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
if (i %in% error_cases) {
monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
std.err = FALSE)
}
else {
monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead")
}
}
optimization may not have succeeded
names(monthly_fits) = names(monthly_max)
for (i in 1:12) {
pdf(paste("../plots/monthly_mle_diag/", sprintf("%02d", i), "_monthly_mle_diag.pdf", sep=""),
width=7, height=7)
par(mfrow=c(2,2))
plot(monthly_fits[[i]])
dev.off()
}
gumbel_qq(data.frame(monthly_max[[9]])[,1], "September")

gumbel_qq(data.frame(monthly_max[[12]])[,1], "December")

Fit R largest order statistics
largest_10 = get_monthly(apply.monthly(prod_ts, function(x) order(x, decreasing=T)[1:4]))
largest_2_fit = lapply(largest_10, function(x) rlarg.fit(x[ , 1:2]))
$conv
[1] 0
$nllh
[1] 381.7301
$mle
[1] 172.4434267 63.1813324 -0.8246329
$se
[1] 10.3518367 6.2427031 0.1133745
$conv
[1] 0
$nllh
[1] 366.2861
$mle
[1] 155.5451810 51.4464515 -0.7267874
$se
[1] 8.3567636 4.6071290 0.1019313
$conv
[1] 0
$nllh
[1] 348.4217
$mle
[1] 195.3919929 40.8441452 -0.8324715
$se
[1] 6.33288212 4.00545299 0.08820541
$conv
[1] 0
$nllh
[1] 375.033
$mle
[1] 162.9292046 57.8831308 -0.7755148
$se
[1] 9.282767 5.432561 0.099956
$conv
[1] 0
$nllh
[1] 380.4256
$mle
[1] 157.2651459 61.7442705 -0.6554637
$se
[1] 11.5888901 5.2748342 0.1508781
$conv
[1] 0
$nllh
[1] 364.4628
$mle
[1] 104.8006071 51.1390400 -0.1253138
$se
[1] 9.153534 4.541184 0.147210
$conv
[1] 0
$nllh
[1] 386.7889
$mle
[1] 128.6412577 66.3424033 -0.4883028
$se
[1] 11.877915 5.260236 0.139483
$conv
[1] 0
$nllh
[1] 380.7543
$mle
[1] 130.1181620 62.4551520 -0.3960438
$se
[1] 10.9040496 4.7896609 0.1237002
$conv
[1] 0
$nllh
[1] 361.8412
$mle
[1] 111.0849233 49.8496920 -0.2990619
$se
[1] 8.2812262 3.8324125 0.1013576
$conv
[1] 0
$nllh
[1] 377.358
$mle
[1] 75.1622227 56.9409902 0.2522703
$se
[1] 10.0789585 7.2986557 0.1984848
$conv
[1] 0
$nllh
[1] 376.1717
$mle
[1] 106.7510273 60.2707332 -0.1192325
$se
[1] 10.4205289 5.3759827 0.1333882
$conv
[1] 0
$nllh
[1] 377.2483
$mle
[1] 126.4665127 59.0621082 -0.2939312
$se
[1] 9.9836743 4.7019303 0.1170011
get_se = function(x, ix) {
if (is.null(x$std.err)) 0
else x$std.err[ix]
}
mle_loc = unlist(lapply(monthly_fits, function(x) x$estimate[1]))
mle_loc_se = unlist(lapply(monthly_fits, get_se, 1))
mle_scale = unlist(lapply(monthly_fits, function(x) x$estimate[2]))
mle_scale_se = unlist(lapply(monthly_fits, get_se, 2))
mle_shape = unlist(lapply(monthly_fits, function(x) x$estimate[3]))
mle_shape_se = unlist(lapply(monthly_fits, get_se, 3))
plot_w_err = function(x, y, se, se.conf_mult = 1,title = NULL) {
max_ix = which.max(y)
min_ix = which.min(y)
plot(x, y,
ylim = c(y[min_ix] - se.conf_mult*se[min_ix], y[max_ix] + se.conf_mult*se[max_ix]),
main = title)
arrows(x,y-se.conf_mult*se,x,y+se.conf_mult*se, code=3, length=0.02, angle = 90)
}
plot_w_err(1:12, mle_loc, mle_loc_se, 1, "Location Parameter as Estimated with Likelihood")

plot_w_err(1:12, mle_scale, mle_scale_se, 1, "Scale Parameter as Estimated with Likelihood")

plot_w_err(1:12, mle_shape, mle_shape_se, 1, "Shape Parameter as Estimated with Likelihood")

** Fit with Bayesian Methods for months separately**
# Fits GEV distribution with bayesian method for given parameters
bayes_fitter = function(x,
init = c(1e3, 1e3, 0.1), # Initial values
mat = diag(c(10000,10000,100)),
psd = c(500,0.1,0.1), # Proposed SDev
nit = 3000, # Nb Iterations
thinning = 50, # Thinning Factor
do_diagn = FALSE, # Bool whether to show diagnostic plots
do_autoreg = FALSE, # Bool whether to show autoreg plots
seed = 42 # Seed
) {
set.seed(seed)
pn = prior.norm(mean=c(0,0,0),cov=mat)
post = posterior(nit, init=init, prior=pn, lh="gev", data=x, psd=psd)
if(do_diagn) {
MCMC<-mcmc(post[1:nit %% thinning == 0, ])
plot(MCMC)
}
if(do_autoreg) {
acf(mcmc(post[1:nit %% thinning == 0, ]))
}
list(posterior = post,
acceptance_rate = attr(mcmc(post),"ar"))
}
# Iteratively fits GEV with bayesian methods, until the fit has
# acceptable acceptance rates (i.e. 0.2 < AR < 0.4). If the AR is too high,
# the proposed SDev for the parameter is multiplied with 1.5. If it's too small,
# the proposed SDev is divided by 2. This is repeated until the acceptance rate
# is good for all parameters, or max_it is reached. Then, a final model is fitted
# with more iterations.
bayes_fitter_param_search = function(x,
psd_init = c(500,0.1,0.1), # Initial proposed SDev
nit_full = 3000, # Nb Iterations for final model
nit_search = 150, # Nb Iterations for param search
do_diagn = FALSE, # Bool whether to show diagnostic plots
do_autoreg = FALSE, # Bool whether to show autoreg plots
max_it = 20,
... # Additional params passed to bayes_fitter
)
{
# Iterate until desired acceptance rate is obtained
cont = TRUE
psd = psd_init
it = 0
while(cont) {
it = it+1
if (it > max_it) {
warning("The ")
}
fit = bayes_fitter(x, psd=psd, nit=nit_search, do_diagn=FALSE,
do_autoreg=FALSE,...)
acc_rates = fit$acceptance_rate[1, 1:3]
too_high = acc_rates > .4
too_low = acc_rates < .2
if (all(!too_high) && all(!too_low)) { # All acceptance rates lie within threshold
cont=FALSE
} else if (it > max_it) { # max_it is reached
cont=FALSE
warning("max_it was reached")
} else { # Correct values which have wrong threshold
psd[too_high] = psd[too_high] * 1.5
psd[too_low] = psd[too_low] / 2
}
}
# Fit final model
bayes_fitter(x, psd=psd, nit=nit_full, do_diagn=do_diagn,
do_autoreg=do_autoreg, ...)
}
monthly_bayes_fit = lapply(monthly_max, bayes_fitter_param_search, do_diagn = TRUE,
do_autoreg = TRUE,
psd = c(500,0.3,0.3), nit_full=30000, nit_search = 30000,
thinning = 300)
























acceptance_rates = lapply(monthly_bayes_fit, function(x) x$acceptance_rate[1, ])
print(acceptance_rates)
$jan
mu sigma xi total
0.31 0.32 0.33 0.32
$feb
mu sigma xi total
0.24 0.30 0.33 0.29
$mar
mu sigma xi total
0.24 0.32 0.20 0.25
$apr
mu sigma xi total
0.23 0.35 0.28 0.29
$may
mu sigma xi total
0.24 0.29 0.21 0.25
$jun
mu sigma xi total
0.24 0.32 0.22 0.26
$jul
mu sigma xi total
0.23 0.20 0.29 0.24
$aug
mu sigma xi total
0.24 0.21 0.25 0.23
$sep
mu sigma xi total
0.23 0.27 0.20 0.23
$oct
mu sigma xi total
0.24 0.38 0.23 0.28
$nov
mu sigma xi total
0.24 0.34 0.34 0.31
$dec
mu sigma xi total
0.34 0.33 0.30 0.33
bayes_params = lapply(monthly_bayes_fit, function(x) apply(x$posterior, 2, median))
bayes_stderr = lapply(monthly_bayes_fit, function(x) apply(x$posterior, 2, sd))
TODO-> R largest fit
PART 2 We perform 2 tests, first for linear depenance, then for step dependance (trend = int(rescaled(trend)>0))
First, we check if the location parameter depends on time using a likelihood ratio test
#monthly_fits = lapply(monthly_max,
# function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
trend = 1:length(as.data.frame(monthly_max[[12]])$V1)
trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
print(i)
fit_const = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
std.err = FALSE)
fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
nsloc = trend,
std.err = FALSE)
ratios[[i]] = fit_const$dev-fit_dependant$dev
}
[1] 1
optimization may not have succeeded
[1] 2
optimization may not have succeeded
[1] 3
[1] 4
optimization may not have succeeded
[1] 5
optimization may not have succeeded
[1] 6
[1] 7
[1] 8
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
optimization may not have succeeded
[1] 12
names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for time independance, \nBonferroni multiple Testing", xlab="Month",ylab="Likelihood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")

Now, let’s check for independance from ENSO
#monthly_fits = lapply(monthly_max,
# function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
# split nino data into months
n = nino34
dim(n)=c(12,length(as.data.frame(monthly_max[[12]])$V1))
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
print(i)
trend = n[i,]
trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
fit_const = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
std.err = FALSE)
fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
nsloc = trend,
std.err = FALSE)
ratios[[i]] = fit_const$dev-fit_dependant$dev
}
[1] 1
optimization may not have succeeded
[1] 2
optimization may not have succeeded
[1] 3
[1] 4
optimization may not have succeeded
[1] 5
[1] 6
optimization may not have succeeded
[1] 7
[1] 8
optimization may not have succeeded
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
[1] 12
names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for independance from ENSO, \nBonferroni multiple testing", xlab="Month",ylab="Likelyhood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")

Now the same for step dependance
#monthly_fits = lapply(monthly_max,
# function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
# split nino data into months
n = nino34
dim(n)=c(12,length(as.data.frame(monthly_max[[12]])$V1))
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
print(i)
trend = n[i,]
trend = as.integer(trend>0)
trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
fit_const = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
std.err = FALSE)
fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
nsloc = trend,
std.err = FALSE)
ratios[[i]] = fit_const$dev-fit_dependant$dev
}
[1] 1
[1] 2
[1] 3
optimization may not have succeeded
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
optimization may not have succeeded
[1] 12
names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for independance from ENSO, \nBonferroni multiple testing", xlab="Month",ylab="Likelyhood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")

#monthly_fits = lapply(monthly_max,
# function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
trend = 1:length(as.data.frame(monthly_max[[12]])$V1)
trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
trend = as.integer(trend>0)
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
print(i)
fit_const = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
std.err = FALSE)
fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
nsloc = trend,
std.err = FALSE)
ratios[[i]] = fit_const$dev-fit_dependant$dev
}
[1] 1
[1] 2
optimization may not have succeeded
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
optimization may not have succeeded
[1] 12
names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for time independance, \nBonferroni multiple Testing", xlab="Month",ylab="Likelihood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")

PART 3 We will now analyse temporal clustering of extremes. For this, we will use the exiplot function from the evd library.
# define a function for getting the extremal indices for each month for a given threshold
monthly_eindex <- function(data, threshold_p, r=0){
ei = list()
for (i in 1:length(data)) {
threshold = quantile(as.data.frame(data[[i]])$V1, threshold_p)
ei[[i]]=exi(as.data.frame(data[[i]])$V1, threshold, r)
}
names(ei) = names(data)
return(ei)
}
ei = monthly_eindex(get_monthly(prod_ts), 0.95)
plot(unlist(ei), main="Extremal Index by Month, 95%-Quantile as Threshold", xlab="Month", ylab="Extremal index")

We observe that the extremal index is 0.25-0.45, we can therefore conclude that we have strong dependance of extremes, with the limiting mean cluster size being roughly from 2 to 4. The clustering has no effect for estimaters based on the (monthly) maximum, but the r largest estimator is influenced by it.
PART 4 First, let’s estimate the 10 year return level using point process approach
monthly_fits_pp = list()
monthly_data = get_monthly(prod_ts)
error_cases = c(5, 9)
month_days = c(31,28,31,30,31,30,31,31,30,31,30,31)
for (i in 1:length(monthly_max)) {
print(i)
threshold = quantile(as.data.frame(monthly_data[[i]])$V1, 0.95)
if (i %in% error_cases) {
monthly_fits_pp[[i]] = fpot(as.data.frame(monthly_data[[i]])$V1,
threshold = threshold,
model="pp",
npp = month_days[i]*8,
cmax = TRUE,
r = 1,
std.err = FALSE,
method = "Nelder-Mead")
}
else {
monthly_fits_pp[[i]] = fpot(as.data.frame(monthly_data[[i]])$V1,
threshold = threshold,
model="pp",
npp = month_days[i]*8,
cmax = TRUE,
r = 1,
method = "Nelder-Mead")
}
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
names(monthly_fits_pp) = names(monthly_data)
for(i in 1:12){
par(mfrow=c(2,2))
plot(monthly_fits_pp[[i]])
}












The fit in february has completely failed, and the others are not very good either
We will still estimate the return levels:
return_level = function(x,period=20){
if (is.list(x)) {
loc = x$estimate[[1]]
scale = x$estimate[[2]]
shape = x$estimate[[3]]
}
if (is.vector(x)) {
loc = x[1]
scale = x[2]
shape = x[3]
}
p = 1/period
level = loc + scale*(((-log(1-p))^-shape-1)/shape)
return(level)
}
return_level_20 = lapply(monthly_fits, return_level) # 20 for testing
return_level_100 = lapply(monthly_fits, return_level, period=100)
return_level_50 = lapply(monthly_fits, return_level, period=50)
plot(unlist(return_level_100),main="100 Year Return level, estimated with point process", xlab="Month",ylab="Return Level")

plot(unlist(return_level_50),main="50 Year Return level, estimated with point process", xlab="Month",ylab="Return Level")

TODO: estimat with mcmc Assuming that we have the posterior densities of the markov chains, call theis function to plot return level histograms
return_level_mcmc = function(posterior,period=20, plot=F){
u = mc.quant(posterior,p=1-1/period,lh="gev")
label_mcmc_rl = sprintf("%s Year Return Level",period)
if(plot) hist(u,nclass=20,prob=T,xlab=label_mcmc_rl, main = "Return Level Histogram")
}
lapply(monthly_bayes_fit, function(x) return_level_mcmc(x$posterior, period=50))
$jan
NULL
$feb
NULL
$mar
NULL
$apr
NULL
$may
NULL
$jun
NULL
$jul
NULL
$aug
NULL
$sep
NULL
$oct
NULL
$nov
NULL
$dec
NULL
lapply(monthly_bayes_fit, function(x) return_level_mcmc(x$posterior, period=100))
$jan
NULL
$feb
NULL
$mar
NULL
$apr
NULL
$may
NULL
$jun
NULL
$jul
NULL
$aug
NULL
$sep
NULL
$oct
NULL
$nov
NULL
$dec
NULL
bayes_50yr_retlvl = lapply(bayes_params, return_level, period=50)
bayes_100yr_retlvl = lapply(bayes_params, return_level, period=100)
** Part 6** We now examine the assymptotic dependance between CAPE and SRH. This can be done by using chi-chibar plots. To eliminate the seasonality, we again consider the monthly data
# plot the chi plot for dependance analysis
data_overall_ts = xts(cbind(cape,srh), order.by = rep(dates, each=8))
data_monthly_bivariate = get_monthly(data_overall_ts)
for (i in 1:length(data_monthly_bivariate)) {
cape.local = as.numeric(data_monthly_bivariate[[i]]$cape)
srh.local = as.numeric(data_monthly_bivariate[[i]]$srh)
dat.cape_srh = cbind(cape.local,srh.local);
par(mfrow=c(1,2))
label_chi = sprintf("Chi plot Month %s",i)
label_chibar = sprintf("Chi Bar plot Month %s",i)
chiplot(dat.cape_srh, main1 = label_chi, which=1);
abline(a=0,b=0,col="red")
chiplot(dat.cape_srh, main2 = label_chibar, which=2);
}
It seems that the confidence intervals for chi contain the 0 horizontal line, we can therefore not conclude that cape and srh are dependant.
Part 7 We will now fit different bivariate models to the monthly maxima and compare them using AIC. We will also estimate the dependance using the logistic link. (if the confidence interval for dep does not contain 1, then we can conclude that there is dependance).
# prepare the data
monthly_max.cape = get_monthly(apply.monthly(data_overall_ts$cape, max))
monthly_max.srh = get_monthly(apply.monthly(data_overall_ts$srh, max))
# fit the different models
zeros = 0*(1:12)
AIC_df = data.frame(zeros,zeros,zeros,zeros,zeros,zeros)
dep.estimate = zeros
dep.sd = zeros
bv_evd_models = c("log","alog","neglog","bilog","ct","negbilog")
colnames(AIC_df) = bv_evd_models
bv_fits = list()
for (i in 1:12) {
print(i)
# fit the data
monthly_max_bivariate = cbind(as.numeric(monthly_max.cape[[i]]),
as.numeric(monthly_max.srh[[i]]))
# fit1 = fbvevd(monthly_max_bivariate,model = "log")
# fit2 = fbvevd(monthly_max_bivariate,model = "alog",std.err = FALSE)
# fit3 = fbvevd(monthly_max_bivariate,model="neglog")
# fit4 = fbvevd(monthly_max_bivariate,model="bilog", std.err = FALSE)
# fit5 = fbvevd(monthly_max_bivariate,model="ct", std.err = FALSE)
# fit6 = fbvevd(monthly_max_bivariate,model="negbilog", std.err = FALSE)
fit1 = fbvevd(monthly_max_bivariate,model = "log", method = "Nelder-Mead")
fit2 = fbvevd(monthly_max_bivariate,model = "alog", method = "Nelder-Mead", std.err = FALSE)
fit3 = fbvevd(monthly_max_bivariate,model="neglog", method = "Nelder-Mead")
fit4 = fbvevd(monthly_max_bivariate,model="bilog", method = "Nelder-Mead", std.err = FALSE)
fit5 = fbvevd(monthly_max_bivariate,model="ct", method = "Nelder-Mead", std.err = FALSE)
fit6 = fbvevd(monthly_max_bivariate,model="negbilog", method = "Nelder-Mead", std.err = FALSE)
bv_fits[[i]] = list(fit1, fit2, fit3, fit4, fit5, fit6)
names(bv_fits[[i]]) = bv_evd_models
# plot the fit for diagnostics
par(mfrow=c(3,2))
plot(fit1)
# keep track of results
dep.estimate[i] = as.numeric(fit1$estimate["dep"])
dep.sd[i] = as.numeric(fit1$std.err["dep"])
aics = c(AIC(fit1), AIC(fit2), AIC(fit3), AIC(fit4), AIC(fit5), AIC(fit6))
AIC_df[i,] = aics
}
names(bv_fits) = month_names
# show AICs for model selection
print(AIC_df)
# plot depencance values for depenance analysis.
par(mfrow=c(1,1))
plot_w_err(1:12, dep.estimate, dep.sd, 1.96,"Monthly Dependance Parameter vs Independance (Red)")
abline(a=1,b=0,col="red")
When plotting the fit diagnostics, the following look good: fit1, fit2, fit4, fit5, fit6
The fit could not be plotted for fit3.
We believe that there is no reason to assume that a different model should be used for different months, so we chose the logistic model. It is one of the best models in every month, and if it isn’t, the difference is only very small.
The plot of the dependance parameter shows that our analysis does coincides with what we saw for the chi-chibar plot, the two series seem to have independant extremes.
# Return levels for model with best AUC for every month
bivar_100yr_retlvl = vector()
bivar_50yr_retlvl = vector()
for (i in 1:12) {
n_sims = 500000
mod_ix = c("neglog"=3) #Always use neglog
#mod_ix = which.min(AIC_df[i, ]) # Index of model with best AIC for given month
estimate = bv_fits[[i]][[mod_ix]]$estimate
model=names(mod_ix)
# Pass appropriate parameters according to model
if (model %in% c("log", "neglog")) {
sim_vals = rbvevd(n=n_sims,
dep=estimate["dep"],
mar1=estimate[c("loc1", "scale1", "shape1")],
mar2=estimate[c("loc2", "scale2", "shape2")],
model=model)
}
else if (model %in% c("bilog", "ct", "negbilog")) {
sim_vals = rbvevd(n=n_sims,
alpha=estimate["alpha"],
beta=estimate["beta"],
mar1=estimate[c("loc1", "scale1", "shape1")],
mar2=estimate[c("loc2", "scale2", "shape2")],
model=model)
}
else if (model %in% c("alog")) {
sim_vals = rbvevd(n=n_sims,
dep=estimate["dep"],
asy=c(estimate["asy1"], estimate["asy2"]),
mar1=estimate[c("loc1", "scale1", "shape1")],
mar2=estimate[c("loc2", "scale2", "shape2")],
model=model)
}
# Obtain simulated PROD value, get return levels
sim_prod = sqrt(sim_vals[ , 1]) * sim_vals[ , 2]
# GET RETURN LEVELS WITH MEAN OF MAX, DOESN'T REALLY MAKE SENSE
# #bivar_100yr_retlvl[i] = mean(apply(matrix(sim_prod, ncol=100), 2, max, na.rm=T),
# na.rm=T)
# #bivar_50yr_retlvl[i] = mean(apply(matrix(sim_prod, ncol=50), 2, max, na.rm=T),
# na.rm=T)
# TAKE N-TH LARGEST ORDER STATISTIC, WITH N = N_SIMS/RET_INTERVAL
#bivar_100yr_retlvl[i] = sort(sim_prod, #partial=n_sims/100,
#decreasing=T)[n_sims/100]
#bivar_50yr_retlvl[i] = sort(sim_prod, #partial=n_sims/50,
#decreasing=T)[n_sims/50]
# USE QUANTILES
bivar_100yr_retlvl[i] = quantile(sim_prod, probs=1-1/100, na.rm = T)
bivar_50yr_retlvl[i] = quantile(sim_prod, probs =1- 1/50, na.rm = T)
}
plot(1:12, return_level_100, col="blue", ylim=c(0, 100000),
main="100 Year Return Levels", pch=19)
points(1:12, bivar_100yr_retlvl, col="red", pch=19)
points(1:12, bayes_100yr_retlvl, col="green", pch=19)
legend("topright", c("Bayesian", "Point Process", "Bivariate"),
fill=c("green","blue","red"))
plot(1:12, return_level_50, col="blue", ylim=c(0, 100000),
main="50 Year Return Levels", pch=19)
points(1:12, bivar_50yr_retlvl, col="red", pch=19)
points(1:12, bayes_50yr_retlvl, col="green", pch=19)
legend("topright", c("Bayesian", "Point Process", "Bivariate"),
fill=c("green","blue","red"))
LS0tCnRpdGxlOiAiVGh1bmRlcnN0cm9tIEFuYWx5c2lzIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCioqSW1wb3J0cyoqCmBgYHtyfQpsaWJyYXJ5KGV2ZCk7IApsaWJyYXJ5KGV2ZGJheWVzKTsKbGlicmFyeShjb2RhKTsKbGlicmFyeShpc21ldik7CmxpYnJhcnkoeHRzKTsKbGlicmFyeShnZ3Bsb3QyKTsKYGBgCgoKKipMb2FkaW5nIHRoZSBEYXRhKioKCmBgYHtyfQpsb2FkKCIuLi9kYXRhL0NBUEVfTWluZGVyX1J5Y2hlbmVyX01hbHNvdC5SRGF0YSIpCmxvYWQoIi4uL2RhdGEvTklOTzM0LlJEYXRhIikKbG9hZCgiLi4vZGF0YS9TUkhfTWluZGVyX1J5Y2hlbmVyX01hbHNvdC5SRGF0YSIpCmBgYAoKKipHZW5lcmF0ZSBQUk9EKioKYGBge3J9CnByb2QgPSBzcXJ0KGNhcGUpKnNyaApgYGAKCgoKKiogQ3JlYXRlIFRpbWUgU2VyaWVzIE9iamVjdHMgKioKYGBge3J9CmRhdGVzIDwtIHNlcS5EYXRlKGFzLkRhdGUoIjE5NzktMS0xIiksIGFzLkRhdGUoIjIwMTUtMTItMzEiKSwgYnk9MSkKZmViMjlpeCA8LSBmb3JtYXQoYXMuRGF0ZShkYXRlcyksICIlbS0lZCIpID09ICIwMi0yOSIKZGF0ZXMgPC0gZGF0ZXNbIWZlYjI5aXhdCgpwcm9kX3RzID0geHRzKHByb2QsIG9yZGVyLmJ5ID0gcmVwKGRhdGVzLCBlYWNoPTgpKQpgYGAKCgoqKlBsb3Qgb2YgVGltZSBTZXJpZXMqKgpgYGB7cn0KdHNfcGxvdCA9IGF1dG9wbG90KHByb2RfdHMpICsgeGxhYigiVGltZSIpICsgeWxhYigiUFJPRCIpCmdnc2F2ZSgiLi4vcGxvdHMvZnVsbF90aW1lX3Nlcmllcy5wZGYiLCBwbG90PXRzX3Bsb3QsIHdpZHRoID0gNiwgaGVpZ2h0ID0gMykKZ2dzYXZlKCIuLi9wbG90cy9mdWxsX3RpbWVfc2VyaWVzLnBuZyIsIHBsb3Q9dHNfcGxvdCwgd2lkdGggPSAzLCBoZWlnaHQgPSAxLjUpCgpwcm9kX2Z1bGxfbW9udGhseV9tYXggPSBhcHBseS5tb250aGx5KHByb2RfdHMsIG1heCkKbWF4X3RzX3Bsb3QgPSBhdXRvcGxvdChwcm9kX2Z1bGxfbW9udGhseV9tYXgpICsKICB4bGFiKCJUaW1lIikgKyB5bGFiKCJQUk9EIikgKwogIGdlb21fbGluZShzaXplPS41KQpnZ3NhdmUoIi4uL3Bsb3RzL21vbnRobHlfbWF4X3Nlcmllcy5wZGYiLCBwbG90PW1heF90c19wbG90LCB3aWR0aCA9IDYsIGhlaWdodCA9IDMpCmdnc2F2ZSgiLi4vcGxvdHMvbW9udGhseV9tYXhfc2VyaWVzLnBuZyIsIHBsb3Q9bWF4X3RzX3Bsb3QsIHdpZHRoID0gMywgaGVpZ2h0ID0gMS41KQoKbmlubzM0X3Bsb3QgPSBnZ3Bsb3QoKSArIGdlb21fbGluZShhZXMoaW5kZXgocHJvZF9mdWxsX21vbnRobHlfbWF4KSwgbmlubzM0KSkgKwogIHhsYWIoIlRpbWUiKSArIHlsYWIoIk5JTk8gMy40IikKZ2dzYXZlKCIuLi9wbG90cy9uaW5vMzRfcGxvdC5wZGYiLCBwbG90PW5pbm8zNF9wbG90LCB3aWR0aCA9IDYsIGhlaWdodCA9IDMpCmdnc2F2ZSgiLi4vcGxvdHMvbmlubzM0X3Bsb3QucG5nIiwgcGxvdD1uaW5vMzRfcGxvdCwgd2lkdGggPSAzLCBoZWlnaHQgPSAxLjUpCgoKYGBgCgoKCioqR3JvdXAgTW9udGhseSoqCmBgYHtyfQptb250aF9uYW1lcyA9IGMoImphbiIsImZlYiIsIm1hciIsImFwciIsIm1heSIsImp1biIsImp1bCIsImF1ZyIsInNlcCIsIm9jdCIsIm5vdiIsImRlYyIpCmdldF9tb250aGx5ID0gZnVuY3Rpb24oeCkgewogIG91dHB1dCA9IGxpc3QoKQogIGxlbiA9IG5yb3coeCkKICAKICAjIEdldCBtb250aCBvZiBlbGVtZW50CiAgbW9udGggPSB0aW1lKHgpCiAgbW9udGggPSBnc3ViKCIuLi4uLSIsICIiLCBtb250aCkKICBtb250aCA9IGdzdWIoIi0uLiIsICIiLCBtb250aCkKICBtb250aGxpc3QgPSB1bmlxdWUobW9udGgpCiAgZm9yIChpIGluIDE6MTIpIHsKICAgIG91dHB1dFtbaV1dID0geFttb250aCA9PSBtb250aGxpc3RbaV0sXQogIH0KICBuYW1lcyhvdXRwdXQpID0gbW9udGhfbmFtZXMKICByZXR1cm4ob3V0cHV0KQp9Cgptb250aGx5X21heCA9IGdldF9tb250aGx5KGFwcGx5Lm1vbnRobHkocHJvZF90cywgbWF4KSkKciA9IDIKcl9tb250aGx5X21heCA9IGdldF9tb250aGx5KGFwcGx5Lm1vbnRobHkocHJvZF90cywgZnVuY3Rpb24oeCkgb3JkZXIoeCwgZGVjcmVhc2luZz1UKVsxOnJdKSkKYGBgCgoqKlBsb3QgTW9udGhseSBUaW1lIFNlcmllcyoqCmBgYHtyfQpmb3IgKGkgaW4gMToxMikgewogIHRtcCA9IGF1dG9wbG90KG1vbnRobHlfbWF4W1tpXV0pICsgeGxhYigiVGltZSIpICsgeWxhYigiUFJPRCIpCiAgZ2dzYXZlKHBhc3RlKCIuLi9wbG90cy9tb250aGx5X21heF9zZXJpZXMvIiwgc3ByaW50ZigiJTAyZCIsIGkpLCAibW9udGhseV9tYXhfc2VyaWVzLnBkZiIsIHNlcD0iIiksIAogICAgICAgICBwbG90PXRtcCwgd2lkdGggPSA2LCBoZWlnaHQgPSAzKQp9CmBgYAoKCmBgYHtyfQojIGdldCB0aGUgbW9udGhseSBtYXhpbWEgb2YgcHJvZAptMSA9IGFzLmRhdGEuZnJhbWUoYXBwbHkubW9udGhseShwcm9kX3RzLCBtYXgpKSRWMTsKIyBwcm9kdWNlIHRoZSBndW1iZWwgcXEgcGxvdApndW1iZWxfcXEgPSBmdW5jdGlvbih4LCB0aXRsZT0iIikgcXFwbG90KHFndW1iZWwoYygxOmxlbmd0aCh4KSkvKGxlbmd0aCh4KSsxKSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtYWluID0gcGFzdGUoIkd1bWJlbGwgUS1RIFBsb3QiLCB0aXRsZSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeGxhYiA9ICJUaGVvcmV0aWNhbCBRdWFudGlsZXMiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5bGFiID0gIlNhbXBsZSBRdWFudGlsZXMiKSA7CgpndW1iZWxfcXEobTEpCgojcXFwbG90KHFndW1iZWwoYygxOmxlbmd0aChtMSkpLyhsZW5ndGgobTEpKzEpKSxtMSxtYWluID0gIkd1bWJlbGwgUS1RIFBsb3QiLHhsYWIgPSAiVGhlb3JldGljYWwgUXVhbnRpbGVzIiwgeWxhYiA9ICJTYW1wbGUgUXVhbnRpbGVzIikgOwpgYGAKVGhlIHFxIHBsb3QgZG9lc24ndCBmaXQgdmVyeSB3ZWxsLCBlc3BlY2lhbGx5IGluIHRoZSBsb3dlciB0YWlsLiBUaGlzIGlzIGxpa2VseSBkdWUKdG8gc2Vhc29uYWwgZGVwZW5kZW5jZS4KCioqRml0dGluZyBHRVYgdG8gdGhlIGVudGlyZSBkYXRhKioKYGBge3J9CiMgZml0IGdldmQgd2l0aCBNTEUgYW5kIHByb2R1Y2UgZGlhZ25vc3RpYyBwbG90cwpmaXRtYXguTUxFPC1mZ2V2KG0xLG1ldGhvZD0iTmVsZGVyLU1lYWQiKQpwYXIobWZyb3c9YygyLDIpKQpmaXRtYXguTUxFCnBsb3QoZml0bWF4Lk1MRSkKYGBgClBvb3IgZml0LCBwcm9iYWJseSBiZWNhdXNlIHRoZSBkaXN0cmlidXRpb24gaXNuJ3Qgc3RhdGlvbmFyeS4gVGhpcyBpcyBlc3BlY2lhbGx5IAp2aXNpYmxlIGluIHRoZSBQcm9iYWJpbGl0eSBwbG90LCBpbiB3aGljaCB0aGUgY29uZmlkZW5jZSBiYW5kIGlzIHN1cnBhc3NlZCwgCmluZGljYXRpbmcgYSBwb29yIGZpdC4KCgpgYGB7cn0KIyBmaXQgZ2V2ZCB3aXRoIEJheWVzaWFuIFRlY2huaXF1ZXMKIyB1c2UgdGhlIHByZXZpb3VzIG91dHB1dHMgKHJvdW5kZWQpIGFzIGluaXRpYWwgdmFsdWVzICh1c2UgYSBkaWZmZXJlbnQgc2hhcGUpCmluaXQ8LWMoMS42ZTQsNGUzLDAuMSkKIyBhcmJpdHJhcnkgcHJpb3JzCm1hdCA8LSBkaWFnKGMoMTAwMDAsMTAwMDAsMTAwKSkgCnBuIDwtIHByaW9yLm5vcm0obWVhbj1jKDAsMCwwKSxjb3Y9bWF0KQojIHByb3Bvc2FsIHN0YW5kYXJkIGRldmlhdGlvbiAoZm91bmQgYnkgdHJpYWwgYW5kIGVycm9yIHRvIGdldCAyMCU8YWNjZXB0YW5jZSByYXRlPDQwJSkKcHNkPC1jKDUwMCwwLjAzLDAuMDIpCiMgZHJhdyAzayBzYW1wbGVzIGZyb20gbWFya292IGNoYWluCm5pdCA9IDMwMDAwCnRoaW5uaW5nID0gMzAwCnBvc3Q8LXBvc3RlcmlvcihuaXQsIGluaXQ9aW5pdCxwcmlvcj1wbixsaD0iZ2V2IixkYXRhPW0xLHBzZD1wc2QpCiMgZGlhZ25vc3RpYyBwbG90cwpNQ01DPC1tY21jKHBvc3RbMTpuaXQgJSUgdGhpbm5pbmcgPT0gMCwgXSwgdGhpbj0zMDApIApwbG90KE1DTUMpIAphdHRyKG1jbWMocG9zdCksImFyIikKCmBgYAoKCmBgYHtyfQojTUNNQ19zdGFiIDwtIG1jbWMocG9zdCwgdGhpbj01MCwgc3RhcnQ9MTAwMCkKYWNmKG1jbWMocG9zdFsxOm5pdCAlJSB0aGlubmluZyA9PSAwLCBdKSkKYGBgCldlIG9ic2VydmUgdGhhdCB0aGVyZSBzZWVtIHRvIGJlIG5vIHN1YnN0YW50aWFsIGlzc3VlcyBmcm9tIHRoZSBhdXRvY29ycmVsYXRpb24gCnBsb3RzLiAKCmBgYHtyfQphcHBseShtY21jKHBvc3RbMTpuaXQgJSUgdGhpbm5pbmcgPT0gMCwgXSksMixtZWFuKQphcHBseShtY21jKHBvc3RbMTpuaXQgJSUgdGhpbm5pbmcgPT0gMCwgXSksMixzZCkKCmBgYAoKKiogRml0IHdpdGggTUxFIGZvciBtb250aHMgc2VwYXJhdGVseSoqCmBgYHtyfQojbW9udGhseV9maXRzID0gbGFwcGx5KG1vbnRobHlfbWF4LCAKIyAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSBmZ2V2KGRhdGEuZnJhbWUoeClbLDFdLCBtZXRob2Q9Ik5lbGRlci1NZWFkIikpCm1vbnRobHlfZml0cyA9IGxpc3QoKQplcnJvcl9jYXNlcyA9IGMoOSwgMTIpCmZvciAoaSBpbiAxOmxlbmd0aChtb250aGx5X21heCkpIHsKICBpZiAoaSAlaW4lIGVycm9yX2Nhc2VzKSB7CiAgICBtb250aGx5X2ZpdHNbW2ldXSA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzdGQuZXJyID0gRkFMU0UpCiAgfQogIGVsc2UgewogICAgbW9udGhseV9maXRzW1tpXV0gPSBmZ2V2KGFzLmRhdGEuZnJhbWUobW9udGhseV9tYXhbW2ldXSkkVjEsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1ldGhvZCA9ICJOZWxkZXItTWVhZCIpCiAgfQp9Cm5hbWVzKG1vbnRobHlfZml0cykgPSBuYW1lcyhtb250aGx5X21heCkKYGBgCgpgYGB7cn0KZm9yIChpIGluIDE6MTIpIHsKICBwZGYocGFzdGUoIi4uL3Bsb3RzL21vbnRobHlfbWxlX2RpYWcvIiwgc3ByaW50ZigiJTAyZCIsIGkpLCAiX21vbnRobHlfbWxlX2RpYWcucGRmIiwgc2VwPSIiKSwKICAgICAgd2lkdGg9NywgaGVpZ2h0PTcpCiAgcGFyKG1mcm93PWMoMiwyKSkKICBwbG90KG1vbnRobHlfZml0c1tbaV1dKQogIGRldi5vZmYoKQp9CmBgYAoKCmBgYHtyfQpndW1iZWxfcXEoZGF0YS5mcmFtZShtb250aGx5X21heFtbOV1dKVssMV0sICJTZXB0ZW1iZXIiKQpndW1iZWxfcXEoZGF0YS5mcmFtZShtb250aGx5X21heFtbMTJdXSlbLDFdLCAiRGVjZW1iZXIiKQpgYGAKCioqRml0IFIgbGFyZ2VzdCBvcmRlciBzdGF0aXN0aWNzKioKYGBge3J9CgpsYXJnZXN0XzEwID0gZ2V0X21vbnRobHkoYXBwbHkubW9udGhseShwcm9kX3RzLCBmdW5jdGlvbih4KSBvcmRlcih4LCBkZWNyZWFzaW5nPVQpWzE6NF0pKQpsYXJnZXN0XzJfZml0ID0gbGFwcGx5KGxhcmdlc3RfMTAsIGZ1bmN0aW9uKHgpIHJsYXJnLmZpdCh4WyAsIDE6Ml0pKQpgYGAKCmBgYHtyfQpnZXRfc2UgPSBmdW5jdGlvbih4LCBpeCkgewogIGlmIChpcy5udWxsKHgkc3RkLmVycikpIDAKICBlbHNlIHgkc3RkLmVycltpeF0KfQptbGVfbG9jID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGZ1bmN0aW9uKHgpIHgkZXN0aW1hdGVbMV0pKQptbGVfbG9jX3NlID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGdldF9zZSwgMSkpCm1sZV9zY2FsZSA9IHVubGlzdChsYXBwbHkobW9udGhseV9maXRzLCBmdW5jdGlvbih4KSB4JGVzdGltYXRlWzJdKSkKbWxlX3NjYWxlX3NlID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGdldF9zZSwgMikpCm1sZV9zaGFwZSA9IHVubGlzdChsYXBwbHkobW9udGhseV9maXRzLCBmdW5jdGlvbih4KSB4JGVzdGltYXRlWzNdKSkKbWxlX3NoYXBlX3NlID0gdW5saXN0KGxhcHBseShtb250aGx5X2ZpdHMsIGdldF9zZSwgMykpCmBgYAoKCmBgYHtyfQpwbG90X3dfZXJyID0gZnVuY3Rpb24oeCwgeSwgc2UsIHNlLmNvbmZfbXVsdCA9IDEsdGl0bGUgPSBOVUxMKSB7CiAgbWF4X2l4ID0gd2hpY2gubWF4KHkpCiAgbWluX2l4ID0gd2hpY2gubWluKHkpCiAgcGxvdCh4LCB5LAogICAgICAgeWxpbSA9IGMoeVttaW5faXhdIC0gc2UuY29uZl9tdWx0KnNlW21pbl9peF0sIHlbbWF4X2l4XSArIHNlLmNvbmZfbXVsdCpzZVttYXhfaXhdKSwKICAgICAgIG1haW4gPSB0aXRsZSkKICBhcnJvd3MoeCx5LXNlLmNvbmZfbXVsdCpzZSx4LHkrc2UuY29uZl9tdWx0KnNlLCBjb2RlPTMsIGxlbmd0aD0wLjAyLCBhbmdsZSA9IDkwKQp9CnBsb3Rfd19lcnIoMToxMiwgbWxlX2xvYywgbWxlX2xvY19zZSwgMSwgIkxvY2F0aW9uIFBhcmFtZXRlciBhcyBFc3RpbWF0ZWQgd2l0aCBMaWtlbGlob29kIikKcGxvdF93X2VycigxOjEyLCBtbGVfc2NhbGUsIG1sZV9zY2FsZV9zZSwgMSwgIlNjYWxlIFBhcmFtZXRlciBhcyBFc3RpbWF0ZWQgd2l0aCBMaWtlbGlob29kIikKcGxvdF93X2VycigxOjEyLCBtbGVfc2hhcGUsIG1sZV9zaGFwZV9zZSwgMSwgIlNoYXBlIFBhcmFtZXRlciBhcyBFc3RpbWF0ZWQgd2l0aCBMaWtlbGlob29kIikKCmBgYAoKKiogRml0IHdpdGggQmF5ZXNpYW4gTWV0aG9kcyBmb3IgbW9udGhzIHNlcGFyYXRlbHkqKgpgYGB7cn0KCgojIEZpdHMgR0VWIGRpc3RyaWJ1dGlvbiB3aXRoIGJheWVzaWFuIG1ldGhvZCBmb3IgZ2l2ZW4gcGFyYW1ldGVycwpiYXllc19maXR0ZXIgPSBmdW5jdGlvbih4LCAKICAgICAgICAgICAgICAgICAgICAgICAgaW5pdCA9IGMoMWUzLCAxZTMsIDAuMSksICMgSW5pdGlhbCB2YWx1ZXMKICAgICAgICAgICAgICAgICAgICAgICAgbWF0ID0gZGlhZyhjKDEwMDAwLDEwMDAwLDEwMCkpLAogICAgICAgICAgICAgICAgICAgICAgICBwc2QgPSBjKDUwMCwwLjEsMC4xKSwgIyBQcm9wb3NlZCBTRGV2CiAgICAgICAgICAgICAgICAgICAgICAgIG5pdCA9IDMwMDAsICMgTmIgSXRlcmF0aW9ucwogICAgICAgICAgICAgICAgICAgICAgICB0aGlubmluZyA9IDUwLCAjIFRoaW5uaW5nIEZhY3RvcgogICAgICAgICAgICAgICAgICAgICAgICBkb19kaWFnbiA9IEZBTFNFLCAjIEJvb2wgd2hldGhlciB0byBzaG93IGRpYWdub3N0aWMgcGxvdHMKICAgICAgICAgICAgICAgICAgICAgICAgZG9fYXV0b3JlZyA9IEZBTFNFLCAjIEJvb2wgd2hldGhlciB0byBzaG93IGF1dG9yZWcgcGxvdHMKICAgICAgICAgICAgICAgICAgICAgICAgc2VlZCA9IDQyICMgU2VlZCAKICAgICAgICAgICAgICAgICAgICAgICAgKSB7CiAgc2V0LnNlZWQoc2VlZCkgICAgICAgICAgICAgICAgCiAgcG4gPSBwcmlvci5ub3JtKG1lYW49YygwLDAsMCksY292PW1hdCkKICBwb3N0ID0gcG9zdGVyaW9yKG5pdCwgaW5pdD1pbml0LCBwcmlvcj1wbiwgbGg9ImdldiIsIGRhdGE9eCwgcHNkPXBzZCkKICAKICBpZihkb19kaWFnbikgewogICAgTUNNQzwtbWNtYyhwb3N0WzE6bml0ICUlIHRoaW5uaW5nID09IDAsIF0pIAogICAgcGxvdChNQ01DKSAKICB9CiAgaWYoZG9fYXV0b3JlZykgewogICAgYWNmKG1jbWMocG9zdFsxOm5pdCAlJSB0aGlubmluZyA9PSAwLCBdKSkKICB9CiAgbGlzdChwb3N0ZXJpb3IgPSBwb3N0LCAKICAgICAgIGFjY2VwdGFuY2VfcmF0ZSA9IGF0dHIobWNtYyhwb3N0KSwiYXIiKSkKfQoKCgojIEl0ZXJhdGl2ZWx5IGZpdHMgR0VWIHdpdGggYmF5ZXNpYW4gbWV0aG9kcywgdW50aWwgdGhlIGZpdCBoYXMgCiMgYWNjZXB0YWJsZSBhY2NlcHRhbmNlIHJhdGVzIChpLmUuIDAuMiA8IEFSIDwgMC40KS4gSWYgdGhlIEFSIGlzIHRvbyBoaWdoLCAKIyB0aGUgcHJvcG9zZWQgU0RldiBmb3IgdGhlIHBhcmFtZXRlciBpcyBtdWx0aXBsaWVkIHdpdGggMS41LiBJZiBpdCdzIHRvbyBzbWFsbCwKIyB0aGUgcHJvcG9zZWQgU0RldiBpcyBkaXZpZGVkIGJ5IDIuIFRoaXMgaXMgcmVwZWF0ZWQgdW50aWwgdGhlIGFjY2VwdGFuY2UgcmF0ZQojIGlzIGdvb2QgZm9yIGFsbCBwYXJhbWV0ZXJzLCBvciBtYXhfaXQgaXMgcmVhY2hlZC4gVGhlbiwgYSBmaW5hbCBtb2RlbCBpcyBmaXR0ZWQKIyB3aXRoIG1vcmUgaXRlcmF0aW9ucy4gCmJheWVzX2ZpdHRlcl9wYXJhbV9zZWFyY2ggPSBmdW5jdGlvbih4LAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcHNkX2luaXQgPSBjKDUwMCwwLjEsMC4xKSwgIyBJbml0aWFsIHByb3Bvc2VkIFNEZXYKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG5pdF9mdWxsID0gMzAwMCwgIyBOYiBJdGVyYXRpb25zIGZvciBmaW5hbCBtb2RlbAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbml0X3NlYXJjaCA9IDE1MCwgIyBOYiBJdGVyYXRpb25zIGZvciBwYXJhbSBzZWFyY2gKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGRvX2RpYWduID0gRkFMU0UsICMgQm9vbCB3aGV0aGVyIHRvIHNob3cgZGlhZ25vc3RpYyBwbG90cwogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZG9fYXV0b3JlZyA9IEZBTFNFLCAjIEJvb2wgd2hldGhlciB0byBzaG93IGF1dG9yZWcgcGxvdHMKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1heF9pdCA9IDIwLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgLi4uICMgQWRkaXRpb25hbCBwYXJhbXMgcGFzc2VkIHRvIGJheWVzX2ZpdHRlcgogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgKSAKewogICMgSXRlcmF0ZSB1bnRpbCBkZXNpcmVkIGFjY2VwdGFuY2UgcmF0ZSBpcyBvYnRhaW5lZAogIGNvbnQgPSBUUlVFCiAgcHNkID0gcHNkX2luaXQKICBpdCA9IDAKICB3aGlsZShjb250KSB7CiAgICBpdCA9IGl0KzEKICAgIGlmIChpdCA+IG1heF9pdCkgewogICAgICB3YXJuaW5nKCJUaGUgIikKICAgIH0KICAgIGZpdCA9IGJheWVzX2ZpdHRlcih4LCBwc2Q9cHNkLCBuaXQ9bml0X3NlYXJjaCwgZG9fZGlhZ249RkFMU0UsIAogICAgICAgICAgICAgICAgICAgICAgIGRvX2F1dG9yZWc9RkFMU0UsLi4uKQogICAgYWNjX3JhdGVzID0gZml0JGFjY2VwdGFuY2VfcmF0ZVsxLCAxOjNdCiAgICAKICAgIHRvb19oaWdoID0gYWNjX3JhdGVzID4gLjQKICAgIHRvb19sb3cgPSBhY2NfcmF0ZXMgPCAuMgogICAgCiAgICBpZiAoYWxsKCF0b29faGlnaCkgJiYgYWxsKCF0b29fbG93KSkgeyAjIEFsbCBhY2NlcHRhbmNlIHJhdGVzIGxpZSB3aXRoaW4gdGhyZXNob2xkCiAgICAgIGNvbnQ9RkFMU0UKICAgIH0gZWxzZSBpZiAoaXQgPiBtYXhfaXQpIHsgIyBtYXhfaXQgaXMgcmVhY2hlZAogICAgICBjb250PUZBTFNFCiAgICAgIHdhcm5pbmcoIm1heF9pdCB3YXMgcmVhY2hlZCIpCiAgICB9IGVsc2UgeyAjIENvcnJlY3QgdmFsdWVzIHdoaWNoIGhhdmUgd3JvbmcgdGhyZXNob2xkCiAgICAgIHBzZFt0b29faGlnaF0gPSBwc2RbdG9vX2hpZ2hdICogMS41CiAgICAgIHBzZFt0b29fbG93XSA9IHBzZFt0b29fbG93XSAvIDIKICAgIH0KICB9CiAgCiAgIyBGaXQgZmluYWwgbW9kZWwKICBiYXllc19maXR0ZXIoeCwgcHNkPXBzZCwgbml0PW5pdF9mdWxsLCBkb19kaWFnbj1kb19kaWFnbiwgCiAgICAgICAgICAgICAgIGRvX2F1dG9yZWc9ZG9fYXV0b3JlZywgLi4uKQogIAp9CiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKCm1vbnRobHlfYmF5ZXNfZml0ID0gbGFwcGx5KG1vbnRobHlfbWF4LCBiYXllc19maXR0ZXJfcGFyYW1fc2VhcmNoLCBkb19kaWFnbiA9IFRSVUUsIAogICAgICAgICAgICAgICAgICAgICAgICAgICBkb19hdXRvcmVnID0gVFJVRSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgcHNkID0gYyg1MDAsMC4zLDAuMyksIG5pdF9mdWxsPTMwMDAwLCBuaXRfc2VhcmNoID0gMzAwMDAsCiAgICAgICAgICAgICAgICAgICAgICAgICAgIHRoaW5uaW5nID0gMzAwKQphY2NlcHRhbmNlX3JhdGVzID0gbGFwcGx5KG1vbnRobHlfYmF5ZXNfZml0LCBmdW5jdGlvbih4KSB4JGFjY2VwdGFuY2VfcmF0ZVsxLCBdKQpwcmludChhY2NlcHRhbmNlX3JhdGVzKQpiYXllc19wYXJhbXMgPSBsYXBwbHkobW9udGhseV9iYXllc19maXQsIGZ1bmN0aW9uKHgpIGFwcGx5KHgkcG9zdGVyaW9yLCAyLCBtZWRpYW4pKQpiYXllc19zdGRlcnIgPSBsYXBwbHkobW9udGhseV9iYXllc19maXQsIGZ1bmN0aW9uKHgpIGFwcGx5KHgkcG9zdGVyaW9yLCAyLCBzZCkpCgoKYGBgCgoKVE9ETy0+IFIgbGFyZ2VzdCBmaXQKCgoKKipQQVJUIDIqKgpXZSBwZXJmb3JtIDIgdGVzdHMsIGZpcnN0IGZvciBsaW5lYXIgZGVwZW5hbmNlLCB0aGVuIGZvciBzdGVwIGRlcGVuZGFuY2UgKHRyZW5kID0gaW50KHJlc2NhbGVkKHRyZW5kKT4wKSkKCkZpcnN0LCB3ZSBjaGVjayBpZiB0aGUgbG9jYXRpb24gcGFyYW1ldGVyIGRlcGVuZHMgb24gdGltZSB1c2luZyBhIGxpa2VsaWhvb2QgcmF0aW8gdGVzdApgYGB7cn0KI21vbnRobHlfZml0cyA9IGxhcHBseShtb250aGx5X21heCwgCiMgICAgICAgICAgICAgICAgICAgICAgZnVuY3Rpb24oeCkgZmdldihkYXRhLmZyYW1lKHgpWywxXSwgbWV0aG9kPSJOZWxkZXItTWVhZCIpKQpyYXRpb3MgPSBsaXN0KCkKdHJlbmQgPSAxOmxlbmd0aChhcy5kYXRhLmZyYW1lKG1vbnRobHlfbWF4W1sxMl1dKSRWMSkKdHJlbmQgPSAodHJlbmQtbWVhbih0cmVuZCkpL3NkKHRyZW5kKSAjIHNjYWxlIGFuZCBjZW50ZXIgY292YXJpYXRlcyBhcyByZWNvbW1lbmRlZAplcnJvcl9jYXNlcyA9IGMoOSwgMTIpCmZvciAoaSBpbiAxOmxlbmd0aChtb250aGx5X21heCkpIHsKICBwcmludChpKQogIAogIGZpdF9jb25zdCA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzdGQuZXJyID0gRkFMU0UpCiAgZml0X2RlcGVuZGFudCA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBuc2xvYyA9IHRyZW5kLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0ZC5lcnIgPSBGQUxTRSkKICAKICByYXRpb3NbW2ldXSA9IGZpdF9jb25zdCRkZXYtZml0X2RlcGVuZGFudCRkZXYgCn0KbmFtZXMocmF0aW9zKSA9IG5hbWVzKG1vbnRobHlfbWF4KQpjaGlfOTVsZXZlbCA9IHFjaGlzcSgxLTAuMDUvMTIsMSkKCnBsb3QodW5saXN0KHJhdGlvcyksbWFpbj0iOTUlIGNvbmZpZGVuY2UgdGVzdCBmb3IgdGltZSBpbmRlcGVuZGFuY2UsIFxuQm9uZmVycm9uaSBtdWx0aXBsZSBUZXN0aW5nIiwgeGxhYj0iTW9udGgiLHlsYWI9Ikxpa2VsaWhvb2QgUmF0aW8gU3RhdGlzdGljIikKYWJsaW5lKGE9Y2hpXzk1bGV2ZWwsYj0wLGNvbD0icmVkIikKCmBgYAoKTm93LCBsZXQncyBjaGVjayBmb3IgaW5kZXBlbmRhbmNlIGZyb20gRU5TTwpgYGB7cn0KI21vbnRobHlfZml0cyA9IGxhcHBseShtb250aGx5X21heCwgCiMgICAgICAgICAgICAgICAgICAgICAgZnVuY3Rpb24oeCkgZmdldihkYXRhLmZyYW1lKHgpWywxXSwgbWV0aG9kPSJOZWxkZXItTWVhZCIpKQpyYXRpb3MgPSBsaXN0KCkKIyBzcGxpdCBuaW5vIGRhdGEgaW50byBtb250aHMKbiA9IG5pbm8zNApkaW0obik9YygxMixsZW5ndGgoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbMTJdXSkkVjEpKQplcnJvcl9jYXNlcyA9IGMoOSwgMTIpCmZvciAoaSBpbiAxOmxlbmd0aChtb250aGx5X21heCkpIHsKICBwcmludChpKQogIHRyZW5kID0gbltpLF0KICAKICB0cmVuZCA9ICh0cmVuZC1tZWFuKHRyZW5kKSkvc2QodHJlbmQpICMgc2NhbGUgYW5kIGNlbnRlciBjb3ZhcmlhdGVzIGFzIHJlY29tbWVuZGVkCiAgZml0X2NvbnN0ID0gZmdldihhcy5kYXRhLmZyYW1lKG1vbnRobHlfbWF4W1tpXV0pJFYxLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtZXRob2QgPSAiTmVsZGVyLU1lYWQiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0ZC5lcnIgPSBGQUxTRSkKICBmaXRfZGVwZW5kYW50ID0gZmdldihhcy5kYXRhLmZyYW1lKG1vbnRobHlfbWF4W1tpXV0pJFYxLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtZXRob2QgPSAiTmVsZGVyLU1lYWQiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIG5zbG9jID0gdHJlbmQsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc3RkLmVyciA9IEZBTFNFKQogIAogIHJhdGlvc1tbaV1dID0gZml0X2NvbnN0JGRldi1maXRfZGVwZW5kYW50JGRldiAKfQpuYW1lcyhyYXRpb3MpID0gbmFtZXMobW9udGhseV9tYXgpCmNoaV85NWxldmVsID0gcWNoaXNxKDEtMC4wNS8xMiwxKQoKcGxvdCh1bmxpc3QocmF0aW9zKSxtYWluPSI5NSUgY29uZmlkZW5jZSB0ZXN0IGZvciBpbmRlcGVuZGFuY2UgZnJvbSBFTlNPLCBcbkJvbmZlcnJvbmkgbXVsdGlwbGUgdGVzdGluZyIsIHhsYWI9Ik1vbnRoIix5bGFiPSJMaWtlbHlob29kIFJhdGlvIFN0YXRpc3RpYyIpCmFibGluZShhPWNoaV85NWxldmVsLGI9MCxjb2w9InJlZCIpCgpgYGAKCk5vdyB0aGUgc2FtZSBmb3Igc3RlcCBkZXBlbmRhbmNlCmBgYHtyfQojbW9udGhseV9maXRzID0gbGFwcGx5KG1vbnRobHlfbWF4LCAKIyAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSBmZ2V2KGRhdGEuZnJhbWUoeClbLDFdLCBtZXRob2Q9Ik5lbGRlci1NZWFkIikpCnJhdGlvcyA9IGxpc3QoKQojIHNwbGl0IG5pbm8gZGF0YSBpbnRvIG1vbnRocwpuID0gbmlubzM0CmRpbShuKT1jKDEyLGxlbmd0aChhcy5kYXRhLmZyYW1lKG1vbnRobHlfbWF4W1sxMl1dKSRWMSkpCmVycm9yX2Nhc2VzID0gYyg5LCAxMikKZm9yIChpIGluIDE6bGVuZ3RoKG1vbnRobHlfbWF4KSkgewogIHByaW50KGkpCiAgdHJlbmQgPSBuW2ksXQogIHRyZW5kID0gYXMuaW50ZWdlcih0cmVuZD4wKQogIAogIHRyZW5kID0gKHRyZW5kLW1lYW4odHJlbmQpKS9zZCh0cmVuZCkgIyBzY2FsZSBhbmQgY2VudGVyIGNvdmFyaWF0ZXMgYXMgcmVjb21tZW5kZWQKICBmaXRfY29uc3QgPSBmZ2V2KGFzLmRhdGEuZnJhbWUobW9udGhseV9tYXhbW2ldXSkkVjEsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1ldGhvZCA9ICJOZWxkZXItTWVhZCIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc3RkLmVyciA9IEZBTFNFKQogIGZpdF9kZXBlbmRhbnQgPSBmZ2V2KGFzLmRhdGEuZnJhbWUobW9udGhseV9tYXhbW2ldXSkkVjEsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1ldGhvZCA9ICJOZWxkZXItTWVhZCIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbnNsb2MgPSB0cmVuZCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzdGQuZXJyID0gRkFMU0UpCiAgCiAgcmF0aW9zW1tpXV0gPSBmaXRfY29uc3QkZGV2LWZpdF9kZXBlbmRhbnQkZGV2IAp9Cm5hbWVzKHJhdGlvcykgPSBuYW1lcyhtb250aGx5X21heCkKY2hpXzk1bGV2ZWwgPSBxY2hpc3EoMS0wLjA1LzEyLDEpCgpwbG90KHVubGlzdChyYXRpb3MpLG1haW49Ijk1JSBjb25maWRlbmNlIHRlc3QgZm9yIGluZGVwZW5kYW5jZSBmcm9tIEVOU08sIFxuQm9uZmVycm9uaSBtdWx0aXBsZSB0ZXN0aW5nIiwgeGxhYj0iTW9udGgiLHlsYWI9Ikxpa2VseWhvb2QgUmF0aW8gU3RhdGlzdGljIikKYWJsaW5lKGE9Y2hpXzk1bGV2ZWwsYj0wLGNvbD0icmVkIikKCmBgYAoKCmBgYHtyfQojbW9udGhseV9maXRzID0gbGFwcGx5KG1vbnRobHlfbWF4LCAKIyAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSBmZ2V2KGRhdGEuZnJhbWUoeClbLDFdLCBtZXRob2Q9Ik5lbGRlci1NZWFkIikpCnJhdGlvcyA9IGxpc3QoKQp0cmVuZCA9IDE6bGVuZ3RoKGFzLmRhdGEuZnJhbWUobW9udGhseV9tYXhbWzEyXV0pJFYxKQp0cmVuZCA9ICh0cmVuZC1tZWFuKHRyZW5kKSkvc2QodHJlbmQpICMgc2NhbGUgYW5kIGNlbnRlciBjb3ZhcmlhdGVzIGFzIHJlY29tbWVuZGVkCnRyZW5kID0gYXMuaW50ZWdlcih0cmVuZD4wKQplcnJvcl9jYXNlcyA9IGMoOSwgMTIpCmZvciAoaSBpbiAxOmxlbmd0aChtb250aGx5X21heCkpIHsKICBwcmludChpKQogIAogIGZpdF9jb25zdCA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzdGQuZXJyID0gRkFMU0UpCiAgZml0X2RlcGVuZGFudCA9IGZnZXYoYXMuZGF0YS5mcmFtZShtb250aGx5X21heFtbaV1dKSRWMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBuc2xvYyA9IHRyZW5kLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0ZC5lcnIgPSBGQUxTRSkKICAKICByYXRpb3NbW2ldXSA9IGZpdF9jb25zdCRkZXYtZml0X2RlcGVuZGFudCRkZXYgCn0KbmFtZXMocmF0aW9zKSA9IG5hbWVzKG1vbnRobHlfbWF4KQpjaGlfOTVsZXZlbCA9IHFjaGlzcSgxLTAuMDUvMTIsMSkKCnBsb3QodW5saXN0KHJhdGlvcyksbWFpbj0iOTUlIGNvbmZpZGVuY2UgdGVzdCBmb3IgdGltZSBpbmRlcGVuZGFuY2UsIFxuQm9uZmVycm9uaSBtdWx0aXBsZSBUZXN0aW5nIiwgeGxhYj0iTW9udGgiLHlsYWI9Ikxpa2VsaWhvb2QgUmF0aW8gU3RhdGlzdGljIikKYWJsaW5lKGE9Y2hpXzk1bGV2ZWwsYj0wLGNvbD0icmVkIikKCmBgYAoKCgoqKlBBUlQgMyoqCldlIHdpbGwgbm93IGFuYWx5c2UgdGVtcG9yYWwgY2x1c3RlcmluZyBvZiBleHRyZW1lcy4gRm9yIHRoaXMsIHdlIHdpbGwgdXNlIHRoZSBleGlwbG90IGZ1bmN0aW9uIGZyb20gdGhlIGV2ZCBsaWJyYXJ5LgoKYGBge3J9CiMgZGVmaW5lIGEgZnVuY3Rpb24gZm9yIGdldHRpbmcgdGhlIGV4dHJlbWFsIGluZGljZXMgZm9yIGVhY2ggbW9udGggZm9yIGEgZ2l2ZW4gdGhyZXNob2xkCm1vbnRobHlfZWluZGV4IDwtIGZ1bmN0aW9uKGRhdGEsIHRocmVzaG9sZF9wLCByPTApewogIGVpID0gbGlzdCgpCiAgZm9yIChpIGluIDE6bGVuZ3RoKGRhdGEpKSB7CiAgICB0aHJlc2hvbGQgPSBxdWFudGlsZShhcy5kYXRhLmZyYW1lKGRhdGFbW2ldXSkkVjEsIHRocmVzaG9sZF9wKQogICAgZWlbW2ldXT1leGkoYXMuZGF0YS5mcmFtZShkYXRhW1tpXV0pJFYxLCB0aHJlc2hvbGQsIHIpCiAgfQogIG5hbWVzKGVpKSA9IG5hbWVzKGRhdGEpCgogIHJldHVybihlaSkKfQoKZWkgPSBtb250aGx5X2VpbmRleChnZXRfbW9udGhseShwcm9kX3RzKSwgMC45NSkKcGxvdCh1bmxpc3QoZWkpLCBtYWluPSJFeHRyZW1hbCBJbmRleCBieSBNb250aCwgOTUlLVF1YW50aWxlIGFzIFRocmVzaG9sZCIsIHhsYWI9Ik1vbnRoIiwgeWxhYj0iRXh0cmVtYWwgaW5kZXgiKQpgYGAKCldlIG9ic2VydmUgdGhhdCB0aGUgZXh0cmVtYWwgaW5kZXggaXMgfjAuMjUtfjAuNDUsIHdlIGNhbiB0aGVyZWZvcmUgY29uY2x1ZGUgdGhhdCB3ZSBoYXZlIHN0cm9uZyBkZXBlbmRhbmNlIG9mIGV4dHJlbWVzLCB3aXRoIHRoZSBsaW1pdGluZyBtZWFuIGNsdXN0ZXIgc2l6ZSBiZWluZyByb3VnaGx5IGZyb20gMiB0byA0LiBUaGUgY2x1c3RlcmluZyBoYXMgbm8gZWZmZWN0IGZvciBlc3RpbWF0ZXJzIGJhc2VkIG9uIHRoZSAobW9udGhseSkgbWF4aW11bSwgYnV0IHRoZSByIGxhcmdlc3QgZXN0aW1hdG9yIGlzIGluZmx1ZW5jZWQgYnkgaXQuCgoKKipQQVJUIDQqKgpGaXJzdCwgbGV0J3MgZXN0aW1hdGUgdGhlIDEwIHllYXIgcmV0dXJuIGxldmVsIHVzaW5nIHBvaW50IHByb2Nlc3MgYXBwcm9hY2gKYGBge3J9Cm1vbnRobHlfZml0c19wcCA9IGxpc3QoKQptb250aGx5X2RhdGEgPSBnZXRfbW9udGhseShwcm9kX3RzKQplcnJvcl9jYXNlcyA9IGMoNSwgOSkKbW9udGhfZGF5cyA9IGMoMzEsMjgsMzEsMzAsMzEsMzAsMzEsMzEsMzAsMzEsMzAsMzEpCmZvciAoaSBpbiAxOmxlbmd0aChtb250aGx5X21heCkpIHsKICBwcmludChpKQogIHRocmVzaG9sZCA9IHF1YW50aWxlKGFzLmRhdGEuZnJhbWUobW9udGhseV9kYXRhW1tpXV0pJFYxLCAwLjk1KQogIAogIGlmIChpICVpbiUgZXJyb3JfY2FzZXMpIHsKICAgIG1vbnRobHlfZml0c19wcFtbaV1dID0gZnBvdChhcy5kYXRhLmZyYW1lKG1vbnRobHlfZGF0YVtbaV1dKSRWMSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0aHJlc2hvbGQgPSB0aHJlc2hvbGQsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbW9kZWw9InBwIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBucHAgPSBtb250aF9kYXlzW2ldKjgsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgY21heCA9IFRSVUUsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgciA9IDEsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc3RkLmVyciA9IEZBTFNFLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1ldGhvZCA9ICJOZWxkZXItTWVhZCIpCiAgfQogIGVsc2UgewogICAgbW9udGhseV9maXRzX3BwW1tpXV0gPSBmcG90KGFzLmRhdGEuZnJhbWUobW9udGhseV9kYXRhW1tpXV0pJFYxLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHRocmVzaG9sZCA9IHRocmVzaG9sZCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtb2RlbD0icHAiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIG5wcCA9IG1vbnRoX2RheXNbaV0qOCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBjbWF4ID0gVFJVRSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICByID0gMSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtZXRob2QgPSAiTmVsZGVyLU1lYWQiKQogIH0KfQpuYW1lcyhtb250aGx5X2ZpdHNfcHApID0gbmFtZXMobW9udGhseV9kYXRhKQpmb3IoaSBpbiAxOjEyKXsKICBwYXIobWZyb3c9YygyLDIpKSAKICBwbG90KG1vbnRobHlfZml0c19wcFtbaV1dKQp9CgpgYGAKVGhlIGZpdCBpbiBmZWJydWFyeSBoYXMgY29tcGxldGVseSBmYWlsZWQsIGFuZCB0aGUgb3RoZXJzIGFyZSBub3QgdmVyeSBnb29kIGVpdGhlcgoKV2Ugd2lsbCBzdGlsbCBlc3RpbWF0ZSB0aGUgcmV0dXJuIGxldmVsczoKYGBge3J9CnJldHVybl9sZXZlbCA9IGZ1bmN0aW9uKHgscGVyaW9kPTIwKXsKICBpZiAoaXMubGlzdCh4KSkgewogICAgbG9jID0geCRlc3RpbWF0ZVtbMV1dCiAgICBzY2FsZSA9IHgkZXN0aW1hdGVbWzJdXQogICAgc2hhcGUgPSB4JGVzdGltYXRlW1szXV0KICB9CiAgaWYgKGlzLnZlY3Rvcih4KSkgewogICAgbG9jID0geFsxXQogICAgc2NhbGUgPSB4WzJdCiAgICBzaGFwZSA9IHhbM10KICB9CiAgcCA9IDEvcGVyaW9kCiAgCiAgbGV2ZWwgPSBsb2MgKyBzY2FsZSooKCgtbG9nKDEtcCkpXi1zaGFwZS0xKS9zaGFwZSkKICByZXR1cm4obGV2ZWw9bGV2ZWwpCn0KcmV0dXJuX2xldmVsXzIwID0gbGFwcGx5KG1vbnRobHlfZml0cywgcmV0dXJuX2xldmVsKSAjIDIwIGZvciB0ZXN0aW5nCnJldHVybl9sZXZlbF8xMDAgPSBsYXBwbHkobW9udGhseV9maXRzLCByZXR1cm5fbGV2ZWwsIHBlcmlvZD0xMDApCnJldHVybl9sZXZlbF81MCA9IGxhcHBseShtb250aGx5X2ZpdHMsIHJldHVybl9sZXZlbCwgcGVyaW9kPTUwKQpwbG90KHVubGlzdChyZXR1cm5fbGV2ZWxfMTAwKSxtYWluPSIxMDAgWWVhciBSZXR1cm4gbGV2ZWwsIGVzdGltYXRlZCB3aXRoIHBvaW50IHByb2Nlc3MiLCB4bGFiPSJNb250aCIseWxhYj0iUmV0dXJuIExldmVsIikKcGxvdCh1bmxpc3QocmV0dXJuX2xldmVsXzUwKSxtYWluPSI1MCBZZWFyIFJldHVybiBsZXZlbCwgZXN0aW1hdGVkIHdpdGggcG9pbnQgcHJvY2VzcyIsIHhsYWI9Ik1vbnRoIix5bGFiPSJSZXR1cm4gTGV2ZWwiKQoKYGBgCgoKVE9ETzogZXN0aW1hdCB3aXRoIG1jbWMKQXNzdW1pbmcgdGhhdCB3ZSBoYXZlIHRoZSBwb3N0ZXJpb3IgZGVuc2l0aWVzIG9mIHRoZSBtYXJrb3YgY2hhaW5zLCBjYWxsIHRoZWlzIGZ1bmN0aW9uIHRvIHBsb3QgcmV0dXJuIGxldmVsIGhpc3RvZ3JhbXMKCmBgYHtyfQpyZXR1cm5fbGV2ZWxfbWNtYyA9IGZ1bmN0aW9uKHBvc3RlcmlvcixwZXJpb2Q9MjAsIHBsb3Q9Ril7CiAgdSA9IG1jLnF1YW50KHBvc3RlcmlvcixwPTEtMS9wZXJpb2QsbGg9ImdldiIpCiAgbGFiZWxfbWNtY19ybCA9IHNwcmludGYoIiVzIFllYXIgUmV0dXJuIExldmVsIixwZXJpb2QpCiAgaWYocGxvdCkgaGlzdCh1LG5jbGFzcz0yMCxwcm9iPVQseGxhYj1sYWJlbF9tY21jX3JsLCBtYWluID0gIlJldHVybiBMZXZlbCBIaXN0b2dyYW0iKQp9CgpsYXBwbHkobW9udGhseV9iYXllc19maXQsIGZ1bmN0aW9uKHgpIHJldHVybl9sZXZlbF9tY21jKHgkcG9zdGVyaW9yLCBwZXJpb2Q9NTApKQpsYXBwbHkobW9udGhseV9iYXllc19maXQsIGZ1bmN0aW9uKHgpIHJldHVybl9sZXZlbF9tY21jKHgkcG9zdGVyaW9yLCBwZXJpb2Q9MTAwKSkKCmJheWVzXzUweXJfcmV0bHZsID0gbGFwcGx5KGJheWVzX3BhcmFtcywgcmV0dXJuX2xldmVsLCBwZXJpb2Q9NTApCmJheWVzXzEwMHlyX3JldGx2bCA9IGxhcHBseShiYXllc19wYXJhbXMsIHJldHVybl9sZXZlbCwgcGVyaW9kPTEwMCkKCgpgYGAKCgoKKiogUGFydCA2KioKV2Ugbm93IGV4YW1pbmUgdGhlIGFzc3ltcHRvdGljIGRlcGVuZGFuY2UgYmV0d2VlbiBDQVBFIGFuZCBTUkguIFRoaXMgY2FuIGJlIGRvbmUgYnkgdXNpbmcgY2hpLWNoaWJhciBwbG90cy4gVG8gZWxpbWluYXRlIHRoZSBzZWFzb25hbGl0eSwgd2UgYWdhaW4gY29uc2lkZXIgdGhlIG1vbnRobHkgZGF0YQpgYGB7cn0KIyBwbG90IHRoZSBjaGkgcGxvdCBmb3IgZGVwZW5kYW5jZSBhbmFseXNpcwpkYXRhX292ZXJhbGxfdHMgPSB4dHMoY2JpbmQoY2FwZSxzcmgpLCBvcmRlci5ieSA9IHJlcChkYXRlcywgZWFjaD04KSkKZGF0YV9tb250aGx5X2JpdmFyaWF0ZSA9IGdldF9tb250aGx5KGRhdGFfb3ZlcmFsbF90cykKCmZvciAoaSBpbiAxOmxlbmd0aChkYXRhX21vbnRobHlfYml2YXJpYXRlKSkgewogIGNhcGUubG9jYWwgPSBhcy5udW1lcmljKGRhdGFfbW9udGhseV9iaXZhcmlhdGVbW2ldXSRjYXBlKQogIHNyaC5sb2NhbCA9IGFzLm51bWVyaWMoZGF0YV9tb250aGx5X2JpdmFyaWF0ZVtbaV1dJHNyaCkKICBkYXQuY2FwZV9zcmggPSBjYmluZChjYXBlLmxvY2FsLHNyaC5sb2NhbCk7CiAgcGFyKG1mcm93PWMoMSwyKSkKICBsYWJlbF9jaGkgPSBzcHJpbnRmKCJDaGkgcGxvdCBNb250aCAlcyIsaSkKICBsYWJlbF9jaGliYXIgPSBzcHJpbnRmKCJDaGkgQmFyIHBsb3QgTW9udGggJXMiLGkpCiAgY2hpcGxvdChkYXQuY2FwZV9zcmgsIG1haW4xID0gbGFiZWxfY2hpLCB3aGljaD0xKTsKICBhYmxpbmUoYT0wLGI9MCxjb2w9InJlZCIpCiAgY2hpcGxvdChkYXQuY2FwZV9zcmgsIG1haW4yID0gbGFiZWxfY2hpYmFyLCB3aGljaD0yKTsKfQpgYGAKCkl0IHNlZW1zIHRoYXQgdGhlIGNvbmZpZGVuY2UgaW50ZXJ2YWxzIGZvciBjaGkgY29udGFpbiB0aGUgMCBob3Jpem9udGFsIGxpbmUsIHdlIGNhbiB0aGVyZWZvcmUgbm90IGNvbmNsdWRlIHRoYXQgY2FwZSBhbmQgc3JoIGFyZSBkZXBlbmRhbnQuCgoKCioqUGFydCA3KioKV2Ugd2lsbCBub3cgZml0IGRpZmZlcmVudCBiaXZhcmlhdGUgbW9kZWxzIHRvIHRoZSBtb250aGx5IG1heGltYSBhbmQgY29tcGFyZSB0aGVtIHVzaW5nIEFJQy4gV2Ugd2lsbCBhbHNvIGVzdGltYXRlIHRoZSBkZXBlbmRhbmNlIHVzaW5nIHRoZSBsb2dpc3RpYyBsaW5rLiAoaWYgdGhlIGNvbmZpZGVuY2UgaW50ZXJ2YWwgZm9yIGRlcCBkb2VzIG5vdCBjb250YWluIDEsIHRoZW4gd2UgY2FuIGNvbmNsdWRlIHRoYXQgdGhlcmUgaXMgZGVwZW5kYW5jZSkuCgpgYGB7cn0KIyBwcmVwYXJlIHRoZSBkYXRhCm1vbnRobHlfbWF4LmNhcGUgPSBnZXRfbW9udGhseShhcHBseS5tb250aGx5KGRhdGFfb3ZlcmFsbF90cyRjYXBlLCBtYXgpKQptb250aGx5X21heC5zcmggPSBnZXRfbW9udGhseShhcHBseS5tb250aGx5KGRhdGFfb3ZlcmFsbF90cyRzcmgsIG1heCkpCgojIGZpdCB0aGUgZGlmZmVyZW50IG1vZGVscwp6ZXJvcyA9IDAqKDE6MTIpCkFJQ19kZiA9IGRhdGEuZnJhbWUoemVyb3MsemVyb3MsemVyb3MsemVyb3MsemVyb3MsemVyb3MpCmRlcC5lc3RpbWF0ZSA9IHplcm9zCmRlcC5zZCA9IHplcm9zCgpidl9ldmRfbW9kZWxzID0gYygibG9nIiwiYWxvZyIsIm5lZ2xvZyIsImJpbG9nIiwiY3QiLCJuZWdiaWxvZyIpCmNvbG5hbWVzKEFJQ19kZikgPSBidl9ldmRfbW9kZWxzCgpidl9maXRzID0gbGlzdCgpCmZvciAoaSBpbiAxOjEyKSB7CiAgcHJpbnQoaSkKICAjIGZpdCB0aGUgZGF0YQogIG1vbnRobHlfbWF4X2JpdmFyaWF0ZSA9IGNiaW5kKGFzLm51bWVyaWMobW9udGhseV9tYXguY2FwZVtbaV1dKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBhcy5udW1lcmljKG1vbnRobHlfbWF4LnNyaFtbaV1dKSkKICAjIGZpdDEgPSBmYnZldmQobW9udGhseV9tYXhfYml2YXJpYXRlLG1vZGVsID0gImxvZyIpCiAgIyBmaXQyID0gZmJ2ZXZkKG1vbnRobHlfbWF4X2JpdmFyaWF0ZSxtb2RlbCA9ICJhbG9nIixzdGQuZXJyID0gRkFMU0UpCiAgIyBmaXQzID0gZmJ2ZXZkKG1vbnRobHlfbWF4X2JpdmFyaWF0ZSxtb2RlbD0ibmVnbG9nIikgCiAgIyBmaXQ0ID0gZmJ2ZXZkKG1vbnRobHlfbWF4X2JpdmFyaWF0ZSxtb2RlbD0iYmlsb2ciLCBzdGQuZXJyID0gRkFMU0UpIAogICMgZml0NSA9IGZidmV2ZChtb250aGx5X21heF9iaXZhcmlhdGUsbW9kZWw9ImN0Iiwgc3RkLmVyciA9IEZBTFNFKSAKICAjIGZpdDYgPSBmYnZldmQobW9udGhseV9tYXhfYml2YXJpYXRlLG1vZGVsPSJuZWdiaWxvZyIsIHN0ZC5lcnIgPSBGQUxTRSkKICAKICBmaXQxID0gZmJ2ZXZkKG1vbnRobHlfbWF4X2JpdmFyaWF0ZSxtb2RlbCA9ICJsb2ciLCBtZXRob2QgPSAiTmVsZGVyLU1lYWQiKQogIGZpdDIgPSBmYnZldmQobW9udGhseV9tYXhfYml2YXJpYXRlLG1vZGVsID0gImFsb2ciLCBtZXRob2QgPSAiTmVsZGVyLU1lYWQiLCBzdGQuZXJyID0gRkFMU0UpCiAgZml0MyA9IGZidmV2ZChtb250aGx5X21heF9iaXZhcmlhdGUsbW9kZWw9Im5lZ2xvZyIsICBtZXRob2QgPSAiTmVsZGVyLU1lYWQiKSAKICBmaXQ0ID0gZmJ2ZXZkKG1vbnRobHlfbWF4X2JpdmFyaWF0ZSxtb2RlbD0iYmlsb2ciLCAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwgc3RkLmVyciA9IEZBTFNFKSAKICBmaXQ1ID0gZmJ2ZXZkKG1vbnRobHlfbWF4X2JpdmFyaWF0ZSxtb2RlbD0iY3QiLCAgbWV0aG9kID0gIk5lbGRlci1NZWFkIiwgc3RkLmVyciA9IEZBTFNFKSAKICBmaXQ2ID0gZmJ2ZXZkKG1vbnRobHlfbWF4X2JpdmFyaWF0ZSxtb2RlbD0ibmVnYmlsb2ciLCBtZXRob2QgPSAiTmVsZGVyLU1lYWQiLCBzdGQuZXJyID0gRkFMU0UpCgogIGJ2X2ZpdHNbW2ldXSA9IGxpc3QoZml0MSwgZml0MiwgZml0MywgZml0NCwgZml0NSwgZml0NikKICBuYW1lcyhidl9maXRzW1tpXV0pID0gYnZfZXZkX21vZGVscwogIAogICMgcGxvdCB0aGUgZml0IGZvciBkaWFnbm9zdGljcwogIHBhcihtZnJvdz1jKDMsMikpIAogIHBsb3QoZml0MSkKICAKICAjIGtlZXAgdHJhY2sgb2YgcmVzdWx0cwogIGRlcC5lc3RpbWF0ZVtpXSA9IGFzLm51bWVyaWMoZml0MSRlc3RpbWF0ZVsiZGVwIl0pCiAgZGVwLnNkW2ldID0gYXMubnVtZXJpYyhmaXQxJHN0ZC5lcnJbImRlcCJdKQogIGFpY3MgPSBjKEFJQyhmaXQxKSwgQUlDKGZpdDIpLCBBSUMoZml0MyksIEFJQyhmaXQ0KSwgQUlDKGZpdDUpLCBBSUMoZml0NikpCiAgQUlDX2RmW2ksXSA9IGFpY3MKCn0KbmFtZXMoYnZfZml0cykgPSBtb250aF9uYW1lcwoKIyBzaG93IEFJQ3MgZm9yIG1vZGVsIHNlbGVjdGlvbgpwcmludChBSUNfZGYpCgojIHBsb3QgZGVwZW5jYW5jZSB2YWx1ZXMgZm9yIGRlcGVuYW5jZSBhbmFseXNpcy4KcGFyKG1mcm93PWMoMSwxKSkKcGxvdF93X2VycigxOjEyLCBkZXAuZXN0aW1hdGUsIGRlcC5zZCwgMS45NiwiTW9udGhseSBEZXBlbmRhbmNlIFBhcmFtZXRlciB2cyBJbmRlcGVuZGFuY2UgKFJlZCkiKQphYmxpbmUoYT0xLGI9MCxjb2w9InJlZCIpCgpgYGAKCldoZW4gcGxvdHRpbmcgdGhlIGZpdCBkaWFnbm9zdGljcywgdGhlIGZvbGxvd2luZyBsb29rIGdvb2Q6CmZpdDEsIGZpdDIsIGZpdDQsIGZpdDUsIGZpdDYKClRoZSBmaXQgY291bGQgbm90IGJlIHBsb3R0ZWQgZm9yIGZpdDMuCgpXZSBiZWxpZXZlIHRoYXQgdGhlcmUgaXMgbm8gcmVhc29uIHRvIGFzc3VtZSB0aGF0IGEgZGlmZmVyZW50IG1vZGVsIHNob3VsZCBiZSB1c2VkIGZvciBkaWZmZXJlbnQgbW9udGhzLCBzbyB3ZSBjaG9zZSB0aGUgbG9naXN0aWMgbW9kZWwuIEl0IGlzIG9uZSBvZiB0aGUgYmVzdCBtb2RlbHMgaW4gZXZlcnkgbW9udGgsIGFuZCBpZiBpdCBpc24ndCwgdGhlIGRpZmZlcmVuY2UgaXMgb25seSB2ZXJ5IHNtYWxsLgoKVGhlIHBsb3Qgb2YgdGhlIGRlcGVuZGFuY2UgcGFyYW1ldGVyIHNob3dzIHRoYXQgb3VyIGFuYWx5c2lzIGRvZXMgY29pbmNpZGVzIHdpdGggd2hhdCB3ZSBzYXcgZm9yIHRoZSBjaGktY2hpYmFyIHBsb3QsIHRoZSB0d28gc2VyaWVzIHNlZW0gdG8gaGF2ZSBpbmRlcGVuZGFudCBleHRyZW1lcy4KCmBgYHtyfQojIFJldHVybiBsZXZlbHMgZm9yIG1vZGVsIHdpdGggYmVzdCBBVUMgZm9yIGV2ZXJ5IG1vbnRoCmJpdmFyXzEwMHlyX3JldGx2bCA9IHZlY3RvcigpCmJpdmFyXzUweXJfcmV0bHZsID0gdmVjdG9yKCkKCmZvciAoaSBpbiAxOjEyKSB7CiAgbl9zaW1zID0gNTAwMDAwCiAgbW9kX2l4ID0gYygibmVnbG9nIj0zKSAjQWx3YXlzIHVzZSBuZWdsb2cKICAjbW9kX2l4ID0gd2hpY2gubWluKEFJQ19kZltpLCBdKSAjIEluZGV4IG9mIG1vZGVsIHdpdGggYmVzdCBBSUMgZm9yIGdpdmVuIG1vbnRoCiAgZXN0aW1hdGUgPSBidl9maXRzW1tpXV1bW21vZF9peF1dJGVzdGltYXRlCiAgbW9kZWw9bmFtZXMobW9kX2l4KQogICMgUGFzcyBhcHByb3ByaWF0ZSBwYXJhbWV0ZXJzIGFjY29yZGluZyB0byBtb2RlbAogIGlmIChtb2RlbCAlaW4lIGMoImxvZyIsICJuZWdsb2ciKSkgewogICAgc2ltX3ZhbHMgPSByYnZldmQobj1uX3NpbXMsCiAgICAgICAgICAgICAgICAgICAgICBkZXA9ZXN0aW1hdGVbImRlcCJdLAogICAgICAgICAgICAgICAgICAgICAgbWFyMT1lc3RpbWF0ZVtjKCJsb2MxIiwgInNjYWxlMSIsICJzaGFwZTEiKV0sCiAgICAgICAgICAgICAgICAgICAgICBtYXIyPWVzdGltYXRlW2MoImxvYzIiLCAic2NhbGUyIiwgInNoYXBlMiIpXSwKICAgICAgICAgICAgICAgICAgICAgIG1vZGVsPW1vZGVsKQoKICB9CiAgZWxzZSBpZiAobW9kZWwgJWluJSBjKCJiaWxvZyIsICJjdCIsICJuZWdiaWxvZyIpKSB7CiAgICBzaW1fdmFscyA9IHJidmV2ZChuPW5fc2ltcywKICAgICAgICAgICAgICAgICAgICAgIGFscGhhPWVzdGltYXRlWyJhbHBoYSJdLAogICAgICAgICAgICAgICAgICAgICAgYmV0YT1lc3RpbWF0ZVsiYmV0YSJdLAogICAgICAgICAgICAgICAgICAgICAgbWFyMT1lc3RpbWF0ZVtjKCJsb2MxIiwgInNjYWxlMSIsICJzaGFwZTEiKV0sCiAgICAgICAgICAgICAgICAgICAgICBtYXIyPWVzdGltYXRlW2MoImxvYzIiLCAic2NhbGUyIiwgInNoYXBlMiIpXSwKICAgICAgICAgICAgICAgICAgICAgIG1vZGVsPW1vZGVsKQoKICB9CiAgZWxzZSBpZiAobW9kZWwgJWluJSBjKCJhbG9nIikpIHsKICAgIHNpbV92YWxzID0gcmJ2ZXZkKG49bl9zaW1zLAogICAgICAgICAgICAgICAgICAgICAgZGVwPWVzdGltYXRlWyJkZXAiXSwKICAgICAgICAgICAgICAgICAgICAgIGFzeT1jKGVzdGltYXRlWyJhc3kxIl0sIGVzdGltYXRlWyJhc3kyIl0pLAogICAgICAgICAgICAgICAgICAgICAgbWFyMT1lc3RpbWF0ZVtjKCJsb2MxIiwgInNjYWxlMSIsICJzaGFwZTEiKV0sCiAgICAgICAgICAgICAgICAgICAgICBtYXIyPWVzdGltYXRlW2MoImxvYzIiLCAic2NhbGUyIiwgInNoYXBlMiIpXSwKICAgICAgICAgICAgICAgICAgICAgIG1vZGVsPW1vZGVsKQogIH0KCiAgIyBPYnRhaW4gc2ltdWxhdGVkIFBST0QgdmFsdWUsIGdldCByZXR1cm4gbGV2ZWxzCiAgc2ltX3Byb2QgPSBzcXJ0KHNpbV92YWxzWyAsIDFdKSAqIHNpbV92YWxzWyAsIDJdCiAgCiAgIyBHRVQgUkVUVVJOIExFVkVMUyBXSVRIIE1FQU4gT0YgTUFYLCBET0VTTidUIFJFQUxMWSBNQUtFIFNFTlNFCiAgIyAjYml2YXJfMTAweXJfcmV0bHZsW2ldID0gbWVhbihhcHBseShtYXRyaXgoc2ltX3Byb2QsIG5jb2w9MTAwKSwgMiwgbWF4LCBuYS5ybT1UKSwKICAjICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbmEucm09VCkKICAjICNiaXZhcl81MHlyX3JldGx2bFtpXSA9IG1lYW4oYXBwbHkobWF0cml4KHNpbV9wcm9kLCBuY29sPTUwKSwgMiwgbWF4LCBuYS5ybT1UKSwKICAjICAgICAgICAgICAgICAgICAgICAgICAgICAgICBuYS5ybT1UKQoKICAjIFRBS0UgTi1USCBMQVJHRVNUIE9SREVSIFNUQVRJU1RJQywgV0lUSCBOID0gTl9TSU1TL1JFVF9JTlRFUlZBTAogICNiaXZhcl8xMDB5cl9yZXRsdmxbaV0gPSBzb3J0KHNpbV9wcm9kLCAjcGFydGlhbD1uX3NpbXMvMTAwLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgI2RlY3JlYXNpbmc9VClbbl9zaW1zLzEwMF0KICAjYml2YXJfNTB5cl9yZXRsdmxbaV0gPSBzb3J0KHNpbV9wcm9kLCAjcGFydGlhbD1uX3NpbXMvNTAsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICNkZWNyZWFzaW5nPVQpW25fc2ltcy81MF0KCiAgIyBVU0UgUVVBTlRJTEVTCiAgYml2YXJfMTAweXJfcmV0bHZsW2ldID0gcXVhbnRpbGUoc2ltX3Byb2QsIHByb2JzPTEtMS8xMDAsIG5hLnJtID0gVCkKICBiaXZhcl81MHlyX3JldGx2bFtpXSA9IHF1YW50aWxlKHNpbV9wcm9kLCBwcm9icyA9MS0gMS81MCwgbmEucm0gPSBUKQoKfQpgYGAKCmBgYHtyfQpwbG90KDE6MTIsIHJldHVybl9sZXZlbF8xMDAsIGNvbD0iYmx1ZSIsIHlsaW09YygwLCAxMDAwMDApLAogICAgIG1haW49IjEwMCBZZWFyIFJldHVybiBMZXZlbHMiLCBwY2g9MTkpCnBvaW50cygxOjEyLCBiaXZhcl8xMDB5cl9yZXRsdmwsIGNvbD0icmVkIiwgcGNoPTE5KQpwb2ludHMoMToxMiwgYmF5ZXNfMTAweXJfcmV0bHZsLCBjb2w9ImdyZWVuIiwgcGNoPTE5KQpsZWdlbmQoInRvcHJpZ2h0IiwgYygiQmF5ZXNpYW4iLCAiUG9pbnQgUHJvY2VzcyIsICJCaXZhcmlhdGUiKSwgCiAgICAgICBmaWxsPWMoImdyZWVuIiwiYmx1ZSIsInJlZCIpKQoKcGxvdCgxOjEyLCByZXR1cm5fbGV2ZWxfNTAsIGNvbD0iYmx1ZSIsIHlsaW09YygwLCAxMDAwMDApLAogICAgIG1haW49IjUwIFllYXIgUmV0dXJuIExldmVscyIsIHBjaD0xOSkKcG9pbnRzKDE6MTIsIGJpdmFyXzUweXJfcmV0bHZsLCBjb2w9InJlZCIsIHBjaD0xOSkKcG9pbnRzKDE6MTIsIGJheWVzXzUweXJfcmV0bHZsLCBjb2w9ImdyZWVuIiwgcGNoPTE5KQpsZWdlbmQoInRvcHJpZ2h0IiwgYygiQmF5ZXNpYW4iLCAiUG9pbnQgUHJvY2VzcyIsICJCaXZhcmlhdGUiKSwgCiAgICAgICBmaWxsPWMoImdyZWVuIiwiYmx1ZSIsInJlZCIpKQoKCmBgYAoKCg==